Data-driven prediction of complex crystal structures of dense lithium

Lithium (Li) is a prototypical simple metal at ambient conditions, but exhibits remarkable changes in structural and electronic properties under compression. There has been intense debate about the structure of dense Li, and recent experiments offered fresh evidence for yet undetermined crystalline phases near the enigmatic melting minimum region in the pressure-temperature phase diagram of Li. Here, we report on an extensive exploration of the energy landscape of Li using an advanced crystal structure search method combined with a machine-learning approach, which greatly expands the scale of structure search, leading to the prediction of four complex Li crystal structures containing up to 192 atoms in the unit cell that are energetically competitive with known Li structures. These findings provide a viable solution to the observed yet unidentified crystalline phases of Li, and showcase the predictive power of the global structure search method for discovering complex crystal structures in conjunction with accurate machine learning potentials.

. Convergence tests on VASP calculations. Convergence tests on total energy in VASP calculations for the cI16 and C2/m structures. Total energies as a function of a kinetic energy cutoff (ENCUT) and b grid spacing of k-point sampling (KSPACING), relative to the results obtained using more stringent parameters (ENCUT = 1500 eV and KSPACING = 0.08Å −1 ).

S II. GENERATION OF THE DEEP POTENTIAL MODELS
The DP models were trained using a highly accurate DFT dataset constructed by the concurrent learning scheme named the Deep-Potential Generator (DP-GEN) [8]. Two DP models were trained. The first model (Li-DP-Hyb1) was initialized from the known experimental phases of Li, and used for accelerating the structure searches. After structure searches, the predicted structures were used to inform the dataset construction of the second model (Li-DP-Hyb2). The Li-DP-Hyb2 model, which contains all information of Li-DP-Hyb1, was used to calculate the Gibbs free energy by thermodynamic integration (TI) through Deep Potential molecular dynamics (DPMD).
A. Protocol for the construction of the dataset DP-GEN [8] is an efficient tool to construct the most compact and adequate dataset for a DP model. Beginning with a preliminary dataset, DP-GEN runs iteratively through the training, exploration and labeling processes.
Preliminary Dataset To initiate the DP-GEN run and guarantee stable sampling in DP-GEN iterations, four initial DP models were firstly trained based on a preliminary dataset, which contains two kinds of configurations: 2. Locally perturbed configurations. The locally optimized configurations were then subject to perturbations, with the mean perturbation distance of atoms set to 0.03Å, and the mean perturbation fraction of cell vectors set to 1.0%. For each structure at each pressure point, 10 perturbations were performed, and included in the preliminary dataset.
Totally, 649 and 1089 configurations were included in the preliminary datasets of the Li-DP-Hyb1 and Li-DP-Hyb2, respectively. The labels (i.e. energies, forces, and virial tensors) of all configurations in the preliminary datasets were obtained by high-accuracy single-point DFT calculations.
Training In each DP-GEN iteration, the first process is the training of four models, using the same training dataset and hyper-parameters, but different random seeds to initialize the model parameters. To enhance the representability of the DP model, we utilized the three-body embedding descriptor and hybridizes it with the original smooth edition DP descriptor. The technical details of these descriptors can be found elsewhere [9]. The sizes of the 3 hidden layers of the two-body embedding net, three-body embedding net, and the fitting net are (20,40,80), (4,8,16) and (240, 240, 240), respectively. The Adam stochastic gradient descent method [10] with default hyper-parameter settings provided by the TensorFlow package [11] was used to train the DP models. During the training, the loss function was: which measures the accuracy of the model according to the differences between the DFT calculated energyÊ k , forceŝ F k iα and the virial tensorΞ k αβ , and those predicted by the DP model. In Eq. (1), B denotes a mini-batch of datasets, and |B| is the batch size. The superscript k denotes the index of the training data in the mini-batch. Each training datum contains a configuration (including the coordinates of atoms, the box basis vectors, and the element types) and its corresponding labels. Prefactors (p ϵ , p f , p ξ ) are a set of hyper-parameters determining the relative importance of the energy, forces, and virial tensor during the training. The prefactors are gradually adjusted according to the learning rate r l (t), which exponentially decays with the training step t: where r 0 l is the learning rate at the beginning of the training, t d denotes the typical timescale of the learning rate decaying and k d denotes the decay rate. The prefactors vary with the learning rate in the following way: where p α (t) is either of the three pre-factors (p ϵ , p f , p ξ ) at training step t. p start α and p limit α are the pre-factors at the beginning and at an infinitely small learning rate, respectively. In practice, a relatively larger force prefactor at the beginning and relatively balanced prefactors at the end of the training can make the best use of the training datasets and achieve relatively good accuracy [12]. Thus, during the DP-GEN iterations, we set p start ϵ = 0.02, p limit ϵ = 1.00, p start f = 1, 000.00, p limit f = 1.00, p start ξ = 0.00, p limit ξ = 0.00, t d = 5000, k d = 0.95, and r l 0 = 1 × 10 −3 . The total number of training steps is 1, 000, 000.
Exploration In each iteration, the second process is the exploration of configuration space using molecular dynamics. The LAMMPS package [13] compiled with the DeePMD-kit [14] support was employed to perform the isothermicisobaric (N P T ) DPMD [15] simulations for the exploration. The simulations run up to 5 ps. The timestep during the exploration was set to 1 fs. The pressure and temperature in the DPMD simulations are listed in Supplementary  Table S1. For the aP160 and oP192 with large unit cells, 64-atom sub-cells are used as the initial configurations of the DPMD simulations. The number of atoms of other structures during exploration is listed in Supplementary Table  S1.
During the DPMD simulation, we evaluated the error of forces by the deviation of force predictions among the four DP models obtained in the training process. If the maximal deviation of atomic forces on a frame in the MD trajectory lies in between an upper bound and a lower bound, the frame is considered as a candidate configuration. At most 150 candidate configurations were labeled by DFT calculations in each DP-GEN iteration. The upper and lower bound were adjusted with respect to the pressure and are listed in Supplementary Table S1. As the number of iterations increases, the deviation of force predictions by independently trained DP models is gradually reduced. The convergence of DP-GEN is reached when the deviations of force predictions for 99.9% of explored configurations are below the lower bound.
Supplementary Table S1. The number of atoms, pressure, temperature, lower and upper bound of deviations of force predictions, and whether included in the dataset of Li-DP-Hyb1 and Li-DP-Hyb2 during the exploration process in the DP-GEN run. Labeling The labels of the candidate configurations, i.e. the energy, force, and virial tensor, were calculated by DFT using VASP [1] code as described in Sec. SI. The candidate configurations and their labels were then appended to the dataset for the model training in the next iteration.

B. Training and refining of the productive model
After the DP-GEN iterations converged, the productive models were trained with longer training steps with the training dataset generated by DP-GEN. To ensure a high training accuracy, the productive models were firstly trained by 9.60 × 10 7 steps with p start e = 0.02, p limit e = 1.00, p start f = 1000.00, p limit f = 1.00, p start ξ = 0.02, p limit ξ = 1.00, t d = 4.80 × 10 5 , and r l 0 = 1.00 × 10 −3 . Then, the productive models were additionally trained by 4.80 × 10 7 steps with p start e = 10.00, p limit e = 10.00, p start f = 1.00, p limit f = 1.00, p start ξ = 0.02 , p limit ξ = 0.50, t d = 2.40 × 10 5 , and r l 0 = 1.00 × 10 −4 . During the training of the productive models, the weight of the locally optimized configurations in the preliminary dataset was increased by ten times.
To further improve the quality of the model, Li-DP-Hyb2 was then subject to refinement training. Refining datasets were constructed for this purpose. Each refining dataset includes two kinds of configurations: 1. Inaccurate ground-state configurations. The quality of the productive model at zero temperature was accessed by reproducing the enthalpy curves of various Li structures. If the error of DP-predicted enthalpy of a structure at a given pressure is larger than 4 meV/atom, or 2 meV/atom near the enthalpy curve minimum in Supplementary  Fig. S3, the prediction is considered inaccurate, and the labels of the corresponding configurations are calculated by DFT and appended to the refining dataset.
2. Inaccurate MD configurations. The quality of the productive model was also accessed by labeled DPMD trajectories. The trajectories are obtained by DPMD simulations in N P T ensemble using the experimental and predicted structures as the initial configurations. All structures were simulated at 50 K, 100 K, 150 K, 200 K, and 250 K under different pressure conditions. The bcc structure was simulated at 0 GPa and 5 GPa; the f cc structure was simulated at 10 GPa, 20 GPa, and 30 GPa; the cI16 structure was simulated at 45 GPa, 60 GPa, and 70 GPa; the aP160 structure was simulated at 45 GPa; the oP192 structure was simulated at 55 GPa; the oP48 structure was simulated at 60 GPa; the oC88, oC40 and tI20 structures were simulated at 70 GPa. In each thermodynamic condition, DPMD runs up to 5 ps with 1 fs timestep. A frame was labeled every 0.5 ps.
Totally, there are 700 labeled frames in all the trajectories. If the testing error of energy at any trajectories is larger than 2.0, 4.0, 6.0, 8.0, and 10.0 meV/atom, at 50, 100, 150, 200, and 250 K, respectively, we considered the prediction inaccurate. The labels of all frames in all DPMD trajectories were then appended to the refining dataset.
Before each refinement training, the refining dataset was merged with the training dataset. The hyper-parameters during the refinement training were r l 0 = 1.00 × 10 −4 , p start e = 10.00, p limit e = 10.00, p start f = 1.00, p limit f = 1.00, p start ξ = 0.02, p limit ξ = 0.50, t d = 8.00 × 10 4 , and the total number of steps for each refinement training is set to 1.60 × 10 7 . Li-DP-Hyb2 was refined repeatedly until none of the predictions on the enthalpies and MD trajectories were considered inaccurate. Totally, four refining datasets, which consist of nearly 3,000 labeled configurations, were merged into the training dataset of Li-DP-Hyb2 during the refinements. There is no refinement training for Li-DP-Hyb1.

C. Reliability of the Li-DP-Hyb2
The reliability of Li-DP-Hyb2 was accessed by two validations: the prediction on enthalpy curves, and testing on an independent dataset obtained from DPMD with wide thermodynamic conditions.
Prediction of the enthalpy curves. The enthalpies relative to cI16 of various Li structures as a function of pressure are shown in Supplementary Fig. S3. The Li-DP-Hyb2 model can accurately reproduce the enthalpy curves for various Li structures. Especially, at the pressure around the minimum of an enthalpy curve, the errors are well below 1.0 meV/atom. Moreover, the Li-DP-Hyb2 can predict the nearly degenerated enthalpies of oC88 and tI20.
Accuracy on an independent testing dataset. The reliability of the Li-DP-Hyb2 model was also validated by the prediction on an independent testing dataset. The independent testing dataset was constructed in a similar way as in the construction of the refining dataset, but obtained from longer DPMD trajectories (10 ps, 50 frames on each trajectory), and covered a wider range of thermodynamic conditions, as listed in Supplementary Table S2. Totally, 7,000 frames were included in the independent testing dataset. The root mean square error (RMSE) of the predicted energy of each phase on their sampling thermodynamic condition is listed in Supplementary Table S2. The RMSE below 50, 100, 150, 200, and 250 K are 0.68, 0.73, 1.02, 1.21, and 1.39 meV/atom, respectively. The comparisons between the DP-predicted and the DFT-calculated energies of various configurations from 45 to 70 GPa in the independent testing dataset are shown in Supplementary Fig. S4.  Supplementary Fig. S4. DP-predicted energy versus DFT-calculated energy of Li structures in the independent testing dataset. The dash-line indicates the ±2 meV/atom error. a 45GPa, b 55GPa, c 60GPa, d 70GPa.

A. aP160 and similar structures
During the structure searches, several geometrically similar structures as the aP160 were found. These structures include P 2 1 /c (8 atoms/cell), P 2 1 (8 atoms/cell), C2/m (8 atoms/cell), Cmca (16 atoms/cell), and P c (80 atoms/cell) structures. The C2/m structure was proposed in a recent computational study [20]. These structures possess smaller unit cells and higher symmetry than the aP160. They are different from the aP160 in the arrangement manners of the prisms or the magnitude of distortion of atomic chains. The C2/m and aP160 structures are shown in Supplementary  Fig. S5(a) and (b), respectively, and all structures are provided in separate files in the POSCAR format.
The enthalpies of these aP160-like structures were calculated using DFT at a high level of accuracy and plotted as a function of pressure relative to the experimental cI16 in Supplementary Fig. S5(c). These structures show nearly degenerated enthalpies to that of the aP160 at 50 and 55 GPa, but become enthalpically less favorable at higher or lower pressure.
We have also calculated the free energy for the recently proposed C2/m structure [20] at 45 GPa as shown in Supplementary Fig. S5(d). Our results show that the aP160 and the C2/m structures are almost energetically degenerate at 0-300 K and both are less favorable in free energy than the cI16. This contradicts the theoretical results reported in Ref. [20], where the C2/m structure is found to be more stable than the cI16 (by less than 1 meV/atom) at temperatures above ∼200 K. This discrepancy can be attributed to the different choices of parameters in the DFT calculations. Here we adopt more tight parameters to ensure better numerical convergence (see Supplementary Sec. I for computational details for DFT calculations).
An interesting phenomenon observed in the aP160 structure is that the atomic chains are gradually distorted with increasing pressure, similar to that observed during the melting of a material ( Supplementary Fig. S6). This phenomenon is responsible for the favorable enthalpy of aP160 and is not observed in other aP160-like structures due to the constraints of symmetry as well as the small unit cells. The other phenomenon observed in aP160, as well as other structures described below, is the wide distribution of distances of the nearest-neighbor contacts ( Supplementary  Fig. S6). At 50 GPa, the Li...Li distances are between 2.01 -2.47Å. This bond alternation is also found in other structures of Li (e.g. the experimental oC40 and oC88 structures) and is a common feature of Li at high pressure as a result of Peierls distortion [21][22][23]. For the aP160 and oP192, the large unit cell makes DFT-based phonon calculation prohibitively expensive, thus the dynamic stability was examined by DPMD simulations using the N P T ensemble with T = 100 K, and P = 45 and 55 GPa, respectively. The 2 × 2 × 2 supercells (1280 atoms for aP160 and 1536 atoms for oP192) were used in DPMD simulations. Supplementary Fig. S7 shows the mean square displacement (MSD) and the MD trajectories during the MD simulations for the two structures. For both structures, the MSDs are observed to be almost constant during the 1 ns simulations, which proves that the systems are dynamically stable in the solids with all atoms near their equilibrium positions.

S IV. FREE ENERGY CALCULATIONS
The calculation of Gibbs free energy of a Li phase takes a two-step process: 1. Calculating the free energy of a given Li phase relative to a reference system, whose free energy is known. It is realized by the Hamiltonian thermodynamic integration (HTI) at a reference thermodynamic condition, (T 0 , P 0 ), where the dynamic stability of the Li phase should be guaranteed.
2. Calculating the free energy difference of the Li phase between the target thermodynamic condition (T 1 , P 1 ) and the reference condition (T 0 , P 0 ). It is realized by the thermodynamic integration (TI) along an isobar (P 0 = P 1 ) in this work.
The principles of HTI and TI are introduced in Ref. [18]. The HTI and TI are implemented with the workflow in Deep Potential thermodynamic integration (DPTI) code, which automatically generates the scripts for DPMD simulations conducted by LAMMPS [13], and manages the pre-processes, post-processes and numeric integrations. For all DPMD simulations, the timestep was set to 1 fs. Before conducting HTI, a N P T DPMD simulation was conducted to obtain the proper simulation box at (P 0 , T 0 ). Then a N V T DPMD simulation is conducted at T 0 to check whether the obtained simulation box is reasonable. The N P T DPMD simulation runs up to 1 ns, and the N V T DPMD simulation runs up to 200 ps.

A. Hamiltonian thermodynamic integration
For a solid phase, the reference system is the Einstein crystal [24] with the spring constant k = 0.045 eV/Å 2 In the Einstein crystal, atoms are non-interacting. Two intermediate systems are introduced to make the integration path as smooth as possible.
The first intermediate system is the modified Einstein crystal which includes a soft-core Lennard-Jones(LJ) interaction [25]. The potential energy U SC is given by which repels atoms at short distances. We set η = 0. MD simulations at each λ run up to 300 ps. The HTI for solid phases was conducted at T 0 = 150 K, except that for oC88 and oP192 structure at 55 GPa, T 0 = 50 K. The numbers of atoms used in the HTI of cI16, oC88, aP160, tI20, oP48, and oP192 structures are 432, 352, 320, 320, 384, and 768, respectively.
For the HTI of Li liquid phase, the reference system is the ideal gas. The constructions of the two intermediate systems are similar to those in the HTI for solid phases. The first intermediate system switches on the soft-core LJ potential on top of the ideal gas, and the second system additionally switches on DP. The parameters in the soft-core potential in Eq.(4) are set to η = 0.5, σ = 1.95, ϵ = 0.01, α = 0.5, and n = 1.0.
The integration path between the ideal gas and the During the HTI of the liquid phase, The number of atoms used in the liquid phase HTI is 432. MD simulations at each λ run up to 600 ps. The starting temperature T 0 for the liquid phase is set to 400 K.